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We present two methods for integrating forced geodesic equations in the Kerr spacetime. The 
methods can accommodate arbitrary forces. As a test case, we compute inspirals caused by a simple 
drag force, mimicking motion in the presence of gas. We verify that both methods give the same 
results for this simple force. We find that drag generally causes eccentricity to increase throughout 
the inspiral. This is a relativistic effect qualitatively opposite to what is seen in gravitational- 
radiation-driven inspirals, and similar to what others have observed in hydrodynamic simulations 
of gaseous binaries. We provide an analytic explanation by deriving the leading order relativistic 
correction to the Newtonian dynamics. If observed, an increasing eccentricity would thus provide 
clear evidence that the inspiral was occurring in a nonvacuum environment. 

Our two methods are especially useful for evolving orbits in the adiabatic regime. Both use the 
method of osculating orbits, in which each point on the orbit is characterized by the parameters of 
the geodesic with the same instantaneous position and velocity. Both methods describe the orbit in 
terms of the geodesic energy, axial angular momentum. Carter constant, azimuthal phase, and two 
angular variables that increase monotonically and are relativistic generalizations of the eccentric 
anomaly. The two methods differ in their treatment of the orbital phases and the representation 
of the force. In the first method, the geodesic phase and phase constant are evolved together as 
a single orbital phase parameter, and the force is expressed in terms of its components on the 
Kinnersley orthonormal tetrad. In the second method, the phase constants of the geodesic motion 
are evolved separately and the force is expressed in terms of its Boyer-Lindquist components. This 
second approach is a direct generalization of earlier work by Pound and Poisson [l[ for planar forces 
in a Schwarzschild background. 



I. INTRODUCTION 

The two-body problem in relativity when one of the 
bodies is much more massive than the other is of great 
interest both theoretically and astrophysically. In this 
limit, the orbit of the smaller body is approximately 
geodesic on short time scales. Deviations from the 
geodesic trajectory arise from the back-reaction on the 
orbit of the spacetime perturbation created by the object, 
but can also arise from external factors such as gravita- 
tional interactions with other bodies, gaseous material in 
the spacetime and so forth. In all these situations, the 
orbit can be described as a geodesic acted on by a per- 
turbing force, which is in general small. In this article, 
we describe techniques for integrating the Kerr geodesic 
equations in the presence of an arbitrary forcing term, 
which can be applied to any of these problems. 

For the back- reaction on the orbit, the perturbing 
force, called the self-force, is of the order of the mass 
ratio fJ./M and it can be computed by a perturbation 
expansion in this small parameter. Computing the lin- 
earized metric perturbation sourced by the compact ob- 
ject and hence the self-force is not an easy task and it 
has taken more than a decade to solve this problem for a 
nonspinning compact object moving in a Schwarzschild 
background The conventional approach treats the 

compact object as a test mass which leads to a diver- 
gence of the field at the position of the particle and this 



must be dealt with using a regularization procedure. The 
extension to Kerr orbits is underway. The techniques de- 
scribed in this paper will be a useful tool in the future 
for constructing trajectories evolving under gravitational 
radiation-reaction. 

The problem of the motion of two bodies with very dif- 
ferent masses is relevant for present and future gravita- 
tional wave detectors. Systems with mass ratios of 1:100 
(intcrmediate-mass-ratio inspirals) could be detected by 
the advanced generation of ground-based detectors that 
are currently under construction . The proposed space- 
based detector LISA Q is expected to detect '--^ 10 — 100 
extreme- mass-ratio inspiral (EMRI) events per year 
These result from the capture of a compact stellar-mass 
object (a white dwarf, neutron star or black hole) by a 
massive black hole (MBH) from a surrounding cusp of 
stars in a galactic nucleus. The captured object gener- 
ates a large number of gravitational wave cycles while it 
is orbiting in the strong field of the MBH, which makes 
these very good sources to use as probes of strong-field 
gravity For both of these classes of source, techniques 
for evolving the orbit under the influence of both grav- 
itational back-reaction and other perturbing forces are 
essential for constructing accurate waveform templates 
and for understanding how external perturbations can 
leave an imprint on the inspiral trajectory 

We present two implementations that can be used to 
integrate geodesic motion in a Kerr background with an 
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external force. We use the method of osculating elements 
extending previous work [l[ for Schwarzschild orbits to 
the Kerr background. The problem of motion under a 
small perturbation is well studied in celestial mechanics 
and is regularly applied to model the motion of satellites 
and small planets. A geodesic in Newtonian mechanics, 
or relativity, is uniquely characterized cither by the three 
components of the particle position vector, r, and the 
three components of the particle velocity, f , at any time 
or by six orbital constants (three orbital constants of the 
motion and three initial phases). There is a one-to-one 
correspondence between the two characterizations. This 
means that any trajectory can be instantaneously identi- 
fied with a geodesic that has the same values of r and r. 
Of course, at two different instances of time, the geodesies 
will differ, but one can smoothly evolve the geodesic pa- 
rameters to reproduce any nongeodesic trajectory. There 
are several approaches to do so and we describe these in 
the next subsection. 



A. Osculating Elements or variation of constants 

As mentioned above we can describe a bound sta- 
ble geodesic by six parameters, which we denote by /. 
In the nonrelativistic case these parameters are simply 
/ = (r, r), while for geodesic motion in Kerr we can take 
/ = {E,L^,Q,^jjo,(l)o,Xo}- Here E is the energy, the 
azimuthal angular momentum, Q is the Carter constant, 
and the remaining phases are defined in Sec. IIIII below. 

At each instant we can therefore identify the true tra- 
jectory with a corresponding geodesic such that r and f 
are the same. This imposes a particular choice of param- 
eters, J, at each instance of time, and the whole trajec- 
tory is thus described by a sequence in the geodesic phase 
space, e.g., I{t) = {E{t),L,{t),L{t),Mt),Mt),Xo{t)}. 
These are referred to as the osculating orbital elements 
at the osculation epoch t jioj . Another name for this ap- 
proach used in the Hamiltonian description is a variation 
of constants. We preserve the form of the equations of 
motion for a geodesic but slowly vary what used to be 
constants of motion in the unperturbed case. There are 
well known techniques for tackling such problems which 
are widely used in Newtonian celestial mechanics and can 
be extended to the relativistic regime. This was demon- 
strated by Pound and Poisson [l| for the trajectory of a 
particle in a Schwarzschild background under the action 
of (post-Newtonian) radiation reaction. 

When we have a perturbed system of the form 

f = fgco+<5f, (1.1) 

we can describe the perturbed trajectory using the os- 
culating elements referred to the orbits of the geodesic 
system f = fgoo. From the chain rule, any one of the 
osculating elements evolves as 

/ = Vr/-f + V„/-r , (1.2) 



in which the subscripts r and v denote derivatives with 
respect to the orbital position and velocity respectively. 
In the absence of the perturbing force, each osculating 
element is constant, so / = V^/ • f -I- V^,/ • fgoo = 0. The 
perturbation equations thus take the rather simple form 

i^Vyl-Sf. (1.3) 

Given an explicit expression for the perturbing force we 
can integrate these equations. 

The osculating element method can be formulated in 
several different ways. There is freedom in the parame- 
terization of the geodesic solution that is used as a basis 
for deriving the osculating element equations, and in the 
basis used to prescribe the force. It is also possible to 
treat the orbital phase constants either as constants of 
the motion that are evolved explicitly or as part of a total 
phase variable which satisfies new equations that depend 
on the perturbation. We will describe two methods for 
treating the Kerr problem: (i) evolution of E,Lz,Q and 
the full orbital phases with the force prescribed with re- 
spect to the Kinnersley orthonormal tetrad; (ii) evolution 
of the orbital constants of motion E, Lz,Q and the initial 
phases, with the force prescribed by its Boyer-Lindquist 
components. 

In the Hamiltonian approach we start with an unper- 
turbed Hamiltonian, Hq and write the equations of mo- 
tion in terms of the constant canonical coordinates and 
momenta, X",Pa (Hamilton- Jacobi approach), which 
are closely related, if not exactly the same, as the six 
constants of motion introduced above, / flTj . If we can 
describe the perturbation as a small addition SH to the 
unperturbed Hamiltonian, then we can describe the equa- 
tions of motion in the same generalized coordinates and 
momenta, which are no longer constants. The derivatives 
of the perturbation SH give the equations for the evolu- 
tion of X^jPa- Quite often those equations are solved 
iteratively starting with an assumption that the orbit is 
unperturbed in the right-hand side (in the 6H). This is 
similar to the adiabatic solution to the osculating element 
equations which we will describe below. The Hamiltonian 
approach (if it can be formulated) would give equations 
equivalent to approach (ii) mentioned above. 

We note that an obvious method of computing inspi- 
rals is to numerically integrate the second-order forced 
geodesic equations directly, taking the fundamental vari- 
ables to be the Boyer-Lindquist coordinates and their 
derivatives with respect to proper time. The key advan- 
tage of the methods discussed in this paper over second- 
order integrations is that they mesh much more natu- 
rally with the adiabatic approximation and more gener- 
ally with two-time-scale approximation techniques [l^ . 
For cxtremc-mass-ratio inspirals driven by radiation re- 
action, the orbital evolution time scale is much longer 
than the orbital time scale for most of the inspiral, un- 
til the orbit becomes close to the innermost stable orbit. 
The adiabatic approximation to the motion gives the mo- 
tion as an expansion in the ratio of the time scales, and 
then there arc various postadiabatic corrections to this. 



3 



Although it is not possible yet to compute numerically 
the full first-order self-force for generic orbits in Kerr, it 
is possible to compute the averaged, dissipative piece of 
this force, which is sufficient to compute leading-order 
adiabatic inspirals [12[. The two-time-scale expansion 
also allows one to go beyond the adiabatic evolution and 
compute the small, rapidly oscillating perturbations to 
the evolution of the orbital variables, as well as the slow 
secular changes to higher order. 

The two-time-scale method cannot be easily applied to 
the second-order, forced geodesic equations, but it can 
be applied to the equations derived in this paper, as we 
discuss in Sees. HIl and IIIII below. In particular, the oscu- 
lating elements method allows us to explicitly estimate 
the orbital average rate of change of the orbital elements. 
This gives us a physical insight into the effect of a per- 
turbing force on the orbit which is otherwise obscured in 
the integration of the second-order equations of motion. 
Estimation of these secular changes also allows us to con- 
struct the adiabatic evolution of the orbit in the regime 
where it is applicable. 



B. Numerical "kludge" waveform 

Another application for the results described in this 
paper is for the construction of numerical kludge wave- 
forms. The numerical kludge waveform for EMRIs is a 
fast and accurate way to compute the long waveforms 
[l3{ that will be needed for EMRI data analysis. These 
are built in a not entirely consistent way, but the ba- 
sic philosophy is to model the underlying trajectory of 
the inspiralling object as accurately as possible in order 
to obtain the best possible phase match between the true 
and approximate waveforms. The approximation is based 
on geodesic motion in the MBH's spacetime, combined 
with a flat spacetime waveform generation expression. 
In the most recent version of the numerical kludge 14l |. 
the instantaneous geodesic orbit was updated by evolving 
the three constants of the motion E,Lz,Q [I5| only. The 
evolution of the constants was obtained by combining 
post-Newtonian results with fits to numerical fluxes ob- 
tained by solving the Teukolsky equation [ij]. However 
this method of evolving the geodesies is not complete, 
as we described above, since we need to evolve the (ini- 
tial) orbital phases together with the orbital constants 
E,Lz,Q- In particular, the natural (and incorrect) way 
to evolve the phase constants, which is to fix them at 
some initial point, leads to significantly different evolu- 
tions in a time or frequency domain implementation of 
the kludge. The desire to resolve this apparent discrep- 
ancy between the two implementations was one of the 
initial motivations for the work described here. This ar- 
ticle outlines the correct way to evolve geodesies under 
the self-force which could be used to further improve the 
numerical kludge waveforms in both time and frequency 
domain descriptions of them. 



C. Main results and the structure of the paper 



In this paper we will give a detailed description of the 
osculating elements approach applied to an arbitrary per- 
turbing force acting on an object undergoing geodesic 
motion in the Kerr spacetime. As an introduction to 
the three dimensional relativistic problem of perturbed 
geodesic motion we will first consider a toy problem in 
Sec. ini We look at the one-dimensional nonlinear oscil- 
lator acted on by an external force. The external force is 
chosen to have two components: a dissipative part and a 
conservative part (which just redefines the energy of the 
system). As we will see later this problem is a very good 
model for the main problem of perturbed motion in the 
Kerr spacetime. We show how two implementations of 
the osculating elements approach work in this simplified 
model and compare the exact evolution with the adia- 
batic approximation. The second of these two implemen- 
tations [in which we evolve the energy and the initial time 
defined as x{to) = 0] allows us to treat the problem an- 
alytically in terms of Jacobi elliptic functions. This one- 
dimensional example allows the reader to understand the 
main approach which we then extend to the problem of 
forced geodesic motion in the Kerr spacetime in Sec. IIIII 
We start that section with an introduction to our nota- 
tion, before describing the osculating elements approach 
using the Kinnersley tetrad and "Hughes" variables (in 
terms of the orbital constants and the total phase vari- 
ables) .We then describe the forced geodesic equations in 
Boyer-Lindquist coordinates, evolving the orbital con- 
stants and the initial conditions, which is a direct exten- 
sion to Kerr of the Schwarzschild results described in [l[ . 
In both cases, we show how wc can explicitly avoid the 
appearance of an apparent divergence in the osculating 
equations of motion at turning points. 



In Sec. IIVI we illustrate our techniques with a prob- 
lem in which the perturbing force is a "gas-drag" force 
proportional to the velocity of the inspiralling compact 
object. This is a toy model for an object inspiralling 
in a gaseous environment around a MBH. We show that 
the different approaches give identical results, and once 
again compare the exact and adiabatic solutions to the 
problem. The influence of the drag force is to drive the 
inspiral of the object, but it also tends to increase the 
eccentricity of the orbit and decrease the orbital inclina- 
tion. Although we primarily use this problem for illustra- 
tive purposes, the increase in eccentricity is an interesting 
result that could have observational consequences. The 
increase in eccentricity is a purely relativistic effect, and 
is to be expected generically, as we discuss in more detail 
in Appendix |D1 in the context of a drag force acting on 
an object in a Schwarzschild background. 

Wc summarize and discuss our findings in the conclud- 
ing section |Vl Some detailed mathematical calculations 
arc included in additional appendixes. 
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II. A SIMPLE MODEL TO ILLUSTRATE 
METHODS USED: THE PERTURBED 
NONLINEAR OSCILLATOR 

In this section we will study in detail the simple model 
of an anharmonic oscillator subject to an external per- 
turbing force, in order to illustrate and explain in a sim- 
ple context the methods that we use for Kerr inspirals in 
subsequent sections of the paper. 

We take the equation of motion for the position x{t) 
of the oscillator to be 



(2.1) 



Here the frequency of the oscillator is chosen to be unity 
for simplicity, /3 > is a parameter governing the size 
of the nonlinear term, Oext is an externally applied per- 
turbing acceleration, which could be a function of both 
the position and the velocity, and e is a small parameter. 
This simple system is similar in some respects to the sys- 
tem of a point particle in orbit about a Kerr black hole 
and subject to the gravitational self-force. The dimen- 
sionless small parameter e in the system (|2.1[) plays the 
role of the mass ratio in the Kerr case, and the external 
acceleration Ccxt is analogous to the self-force. 



A. Analysis using simple phase and energy 
coordinates on phase space 

Consider initially the situation where the is no external 
acceleration. It is useful for some purposes to use a set of 
phase space coordinates for the nonlinear oscillator which 
eliminate the turning points. We define coordinates a and 
"0, functions of x and v = x, by the equations 



1 2 

-a 
2 



X 

sgn(i) 



1.2 



acos-i/) , 
— sgn(sin?/') 



(2.2a) 

(2.2b) 
(2.2c) 



The expression on the right-hand side of Eq. (|2.2ap is 
just the conserved energy of the system, and a is the 
conserved amplitude of the oscillation. The variable if) 
increases monotonically (but not linearly) with time. The 
equations of motion in these variables are 

d = , (2.3a) 
^ = Vl + /?a2(l + cos2V')/2 . (2.3b) 

Now consider turning on the external force. Then the 
right-hand sides of the equations of motion (|2.3p will ac- 
quire terms proportional to e. If we differentiate the def- 
inition (|2.2b[) of ip with respect to t, insert the result into 
the definition ()2.2ap of a, and solve for ip using also Eq. 
(|2.2cp we obtain 



which explicitly shows the extra forcing term. However 
this term contains an apparent divergence aX ip = Q. The 
divergence is only apparent, since d will be constrained 
to vanish when = 0, because the rate at which the force 
does work will vanish when the velocity of the particle is 
zero. 

To see this explicitly, we substitute the definition 
()2.2bj) of -0 into the definition p.2a|) of a, and solve for x 
to get 



sin + /^a2(l -hcos2 V;)/2 



(2.5) 



Next, we differentiate both sides of Eq. (j2.2ap with re- 
spect to t, and simplify the right-hand side using the 
equation of motion ()2.ip . This gives 



(a -I- (ia^)a = eioext 



(2.6) 



Now using the result (|2.5p for x and substituting into Eq. 
(|2.4p gives the final results 



= ^1 + /^a2(l -Hcos2^)/2 



1 - e 



cos ■i/'Oe 



E^l + /3a^(l + cos2 V')/2 



a(l + ^a2) 
sin Qcxt 



(2.7a) 
(2.7b) 



i/i = - cot -0 + yTT^5a2(TTcos20)72 



(2.4) 



where ac^t{x,v) is evaluated at x ~ acos0, and v ~ 
v{a,ip) given by the expression (|2.5p . 

The final result p.7p now casts the system of differen- 
tial equations entirely in terms of the variables a and ip, 
and as expected there are no divergences. Note however 
that Eq. (|2.4p would show a divergence if one used an 
approximate, orbit-averaged version of d instead of the 
exact expression for d. 

In the analogous problem in Kerr, it is very straight- 
forward to compute the analog of the equation of motion 
(|2.4p which contains the apparent divergence. For numer- 
ical work, this form of the equation would be problem- 
atic, since the right-hand side evaluates to 0/0 at turning 
points. Our goal was to attempt to reformulate the equa- 
tions in Kerr analytically, to achieve a form analogous to 
Eq. (|2.7p . where all the divergences have been removed. 
Although it was not clear a priori that this would be 
possible (because of the complexity of the Kerr-orbit dy- 
namical system), we were successful in finding an explic- 
itly finite form of the equations of motion in both sets of 
variables. 

For the problem that we are really interested in, per- 
turbed geodesies in the Kerr spacetime, it will be espe- 
cially useful to consider the adiabatic limit e of small 
external perturbations. So we consider adiabatic per- 
turbations in the context of our example problem. The 
equations of motion (|2.7p for ip and a can be written in 
the general form 

ijj = w(i/;,a)-}-e.g(i)(V',a) + 0(e2) , (2.8a) 
d = eG^'^\tp,a) + 0{e^) . (2.8b) 
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Here on the right-hand side, all the functions are periodic 
functions of "0 with period 2tt. In Appendix iBl we derive 
the limiting; form of the solutions in the limit e — >■ 0; see 
also Ref. [l^l- The leading order or adiabatic solutions 
are given by the following set of steps: 

1. We define the averaging operation, for any function 
f{ip) of V', by 



{f)a 



Jo "^c^(V.,a) 



(2.9) 



The subscript a on the left hand side is a reminder 
that the averaging operation depends on the value 
of a. 

2. We define the averaged functions 

tu{a) = {Lu{iP,a))a , (2.10) 

and 

GW(a)^(G(i)(^,a)), . (2.11) 

3. We solve a pair of ordinary differential equations in 
the slow time parameter 



t = et 



(2.12) 



for two auxiliary functions x*'°''(f) and a*-°-'(f). This 
pair of ordinary differential equations is 



di 
di 



G(i)(a(")(t)) 



(2.13a) 
(2.13b) 



Note that for this step, one does not need to specify 
a value of e. 

4. We can then write down the adiabatic solutions: 



a{t,e) = a'^°\et) , 



^(i,e) 



X(°)(rf),«^"^(et) 



(0)/ 



(2.14a) 
(2.14b) 



where the function .^(x, i) is defined implicitly by 
the equation 



2L 
2n 



/■2ir di/i 

Jo 



(2.15) 



and satisfies S(x + 27r, a) = S(x, a) + 27r. (The in- 
verse of the mapping S essentially maps the given 
phase space coordinates onto action-angle vari- 
ables.) 



Note that there is an asymmetry in how the forcing 
terms g*-^^ and G*-^^ in Eq. (|2.8p enter into the adiabatic 
solution (|2.14p . The function G'^-*, which drives the en- 
ergy evolution, does enter, but the function g*-^-*, which 
drives the phase evolution, does not enter at all. It infiu- 
ences only the post-l-adiabatic solutions. 

Note also that one cannot obtain the adiabatic solu- 
tions by any simple modification of the original differen- 
tial equations. 



B. Analysis exploiting analytic solution to 
un-forced motion 

It is also possible to find an analytic solution to the un- 
perturbed anharmonic oscillator in terms of elliptic func- 
tions. Equation (|2.2ap can be rearranged to give 



2 



X ~ — [Xj^ ~ X ) [X — X_ 



where we have defined the turning points 



4 = 4 (-l±^Jl + 2E|3 



/3 



(2.16) 



(2.17) 



in terms of the energy E, which is set to be twice the 
conserved quantity on the right-hand side of Eq. ()2.2ap , 
and is related to the amplitude of motion a and the non- 
linearity parameter /3 by 



E^ci' + ^Pa" . (2.18) 

For /3 > 0, all of the solutions are bound and oscillate 
periodically in the interval ~x\ < x < x\. Without 
loss of generality we can set a;(io) = 0, in which case 
Eq. (|2.16p can be rearranged and integrated to give 

dy 1 



^{y^-x^_){x\-y^) ^ 



X F 



sm 












X2 


-xl. 



x\ — x'^_ 



x+ 



±\/|(i-io) 



(2.19) 



Here F((/); fc) denotes the Jacobi elliptic integral of the 
first kind Il6i 



F(0;A:) 



sin (f) 



dx 



v/(l-2;2)(l-Pa;2) 
da 



Vl - /s2 sin^ a 



(2.20) 



In the following, we will denote all elliptic integrals by 
bold capital letters. The inverse of this elliptic integral 
is given by the elliptic function sn(u; k) such that 



F(0;fc) 



cj){u; k) 



[sn(u; k)] 



(2.21) 
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For the solutions (|2.19l) . the parameter k and argument 
u are given by 



yi + 2EP-1 

2V1 + 2EI3 



{l + 2Epf'\t-h) , 



(2.22) 
(2.23) 



The relation between x(t) and the elliptic function is then 



= fcsn(M; k) 



(2.24) 



Solving this gives an expression for that is somewhat 
unsatisfying since using it requires manually flipping the 
sign of X. We can instead get a simpler expression if we 
introduce the additional elliptic function 



dn(it; k) = \/l — k'^sii^{u; k) . 



The result is 



(2.25) 



(2.26) 



where sd(u; k) = sn(u; fc)/dn(M; k). 

We will now derive the osculating element equations 
for the variables u and k. The physical variables E and 
to can be obtained from these simply via 



E = 



2fc2(l - fc2) 



/?(l-2fc2)2 ' 

dE _ 4fc dk 
H ^ ^(l-2fc2)3 ' 



ta 
dt 



2fci 



dfc 



dt VI - 2fc2 dt 



(2.27) 

(2.28) 

(2.29) 
(2.30) 



To derive the equations of motion in the osculating cle- 
ment form we need to differentiate sd(M; k) with respect 
to u and k. This gives 



dsd . , , cn(u; fc) 

— — \u\k) — — T 

du^ ' ' dn2(w;fc) 



(2.31) 



9sd, u cn(M; fc) E[0(u; /c); fc]cn(u; A:) 

Sfc^""' ^ ~ kdn{u;k) ^ fc(l-fc2)dn2(u;fc) 



fc sn(M; fc) 
(1 - fc2)dn(w;fc) 



(2.32) 



where we have introduced the elliptic function cn(u;fc), 
which is defined by the analogue of Eq. (j2.2ip but with 
sin^^ replaced by cos~^, and where E(0; k) is the elliptic 
integral of the second kind [l6j : 



E(0;fc) 



"'"-^ , /I -fc2a;2 
r0 



[ da VY 
Jo 



A;2 sin^ a 



(2.33) 



Since the parameter k depends only on the energy, the 
evolution equation can be derived directly from the equa- 
tion for the energy evolution, which follows by differen- 
tiation of Eq. (|2.2ap and use of Eq. (P?T|) : 



dE 

— = 22;eacxt 
dt 



(2.34) 



The evolution equation for u follows from differentiating 
the orbit equation with respect to time and setting this 
equal to the velocity of the unperturbed orbit, which is 
given by Eq. ^Ml- 




2(l-fc2)fc2 
/3(1 - 2P)2 'du 



(2.35) 



Putting these elements together we find the equations for 
the osculating evolution of the orbit 



eaext(a;,i)(l-2A:2)2 



dk 
dt 

du 

dt " VI - 2A:2 



2 ou 



(2.36) 



1 



eaext(x,i)(l-2fc^)W^(l-fc2) 



1 - 2fc2 + 2fc-^ dsd 
/fc(l-fc2)(i-2F)^^"' ^ + 9fc . 



(2.37) 



where the perturbing force is to be evaluated for the 
geodesic position and velocity. 



^ 2fc2(l-fc2) 



/2fc2(l-fc2)asd 

;3(1 - 2fc2)2 9u 



(2.38) 
(2.39) 



We can now derive the adiabatic approximation to the 
solution of Eqs. (|2.36p and (|2.37p following the steps 
described at the end of Sec. Ill Al Eqs. (|2.36p and 
(|2.37p have the same general form as dtp/dt and da/dt 
in Sec. Ill Al That is, we can write them as 

ii = uj{u,k) + eg'^^\u,k) + 0{e^) (2.40a) 
k = eG^^\u,k) + 0{e^) , (2.40b) 

where we have now redefined the functions w, g^^\ and 
G^^\ By comparing against the formula for m, we find 

w(u,fc) = a;(fc) = (1 - 2fc2)-i/2 ^ (2.41) 
As a result; the averaging operation is greatly simplified: 



{f{n,k))i 



U{k) J, 



du f{u, k) 



(2.42) 



where U{k) is the period in u for the general solution 



U{k) =4F(7r/2,fc) = 4K(fc) . 



(2.43) 
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Here K(fc) is the complete elliptic integral of the first 
kind. Note however that this period depends on fc, 
whereas, in the previous parametrization, the period in 
tp was simply 27r. 

As before, the two functions we wish to average are lo 
and G^^\ Since uj is independent of u, 



UJ = {UJj = OJ 



(2.44) 



To make further progress wc must specify the perturbing 
force. Wc take this to be 



(2.45) 



By substituting the Eqs. ((238| and ([239)) for x(u, k) and 
i(u, k) into this expression, we find, omitting the argu- 
ments for the elliptic functions, all of which depend on 
both u and fc, 



G^'Hu, k) = -7fc(l - fc2)(i _ 2fc2) 



+5k\l - fc2)3/2(i „ 2/fc2)y|sd2^ . (2.46) 

The second term in this expression is symmetric about 
zero, and has period U, so it vanishes under the averag- 
ing operation. The first term is also periodic, but it does 
not vanish under averaging since it is always positive. Re- 
calling the relation between dsd/du and the other elliptic 
functions (|2.3ip . 



the interval {0,U}, so 
du dn'" 



(m + l)(l-fc2) 
(m-^2)(2- fc2) / du dn 





m+2 



(m + 3) / du dn 

"'0 



m+4 



(2.52) 



which after using Eq. ()2.51|) gives 
'2 2-A:2 1 



1 - /fc^l 



(2.53) 

The average can be written as an elliptic integral using 



m 



j du dn = E[0(u; k); k] , (2.54) 



where E(u; k) is the elliptic integral of the second kind 
given in Eq. p.33p . and the amplitude function 4>{u] k) is 
given by Eq. (|2.2ip . This leaves us with the final result 



G^i) ==7(l-2/c2) 



22 -k"^ l\ E(fc) l-fc2 
2,~k 



k J K(fc) 3fc 

(2.55) 

where we have used U = 4K(fc), and we use E(fc) = 
E(7r/2,fc) to denote the complete elliptic integral of the 
second kind. 



and exploiting the following identities, 



9 9 

sn~ + cn 
dn^ + k'^sn^ 

we can rewrite G'^^ as 



(2.47) 

(2.48) 
(2.49) 



G(^) = 7(l-2/c2) 



[(1 - k^) (dn-), - (dn-), 



k - ' \ /k 

(2.50) 

The averaging operations can be reduced to just one in- 
tegral by using the identity [l^ 



du dn™ = 



1 



fc2dn™+i 



{m + l){l - fc2) 
+ (m + 2)(2-fc2) J du dn™- 

— (m 4-3) du dn™^^ 



sn cn 



(2.51) 



The first term in square brackets vanishes on the ends of 



C. Example force 

We will illustrate the techniques described above for an 
oscillator subject to the forcing term given in Eq. ()2.45p . 



Ooxt 



-7X - 



(2.56) 



with e = 10"^, l3 = 0.1, 7 = 0.15, S = 0.2 and initial con- 
ditions a;(0) = 1.0, i(0) — 0. The analytic solutions to 
the un-forced motion, as described in Sees. Ill B] and HT^ 
were found to be essentially identical over the full integra- 
tion time (as we would expect since these are both exact 
solutions to the forced motion). In Fig. [TJ we show the 
significant disagreements between (i) using the analytic 
solution to the un-forced motion as described in Sec. Ill Bl 
(labeled "exact"); and (ii) using the adiabatic approx- 
imation to the evolution, given by Eqs. (|2.40p . (j2.44p 
and (|2.55p (labeled "adiabatic"). 

There was no significant disagreement when it came to 
predicting the scale parameter fc^, but there was signifi- 
cant disagreement in what those formalisms predicted for 
the position x{t), the phase u, and the time-offset to. The 
top panel in Fig. [1] shows that the adiabatic and exact 
positions go completely out of phase around t = 2000, af- 
ter which they then continue to pass in and out of phase 
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2000 4000 6000 8000 10000 



t 




2000 4000 6000 8000 10000 



FIG. 1: Comparison between the exact and adiabatic ap- 
proaches to evolving the orbit, both from Sec. IIIBI The top 
panel shows the exact solution x{t) (solid lines), as well as 
the difference between the exact and adiabatic x{t) (dashed 
lines). The insets show close-up views of the first/last 100s of 
the same curves. The bottom panel shows the disagreement 
between the exact and adiabatic predictions for the phase u 
(solid line) and time-offset to (dashed line). 



with each other. This is to be expected, since the adi- 
abatic solution is only an approximation. The bottom 
panel shows disagreements in both the phase u and the 
time-offset to . These grow to several cycles by the end of 
the integration, while the error in (not shown) remains 
small throughout the integration. This is also to be ex- 
pected. Because of the terms we omit, we would expect 
the error in fc^ to scale like e^, while the error in phase u, 
and correspondingly the error in to from Eq. (j2.29[) . will 
scale like e. 

Given that the exact solutions obtained via the an- 
alytic and phase solutions are the same, the choice of 
which parametrization to use must be made on the ba- 
sis of practicality. The integration of the analytic form 
of the equations is more computationally expensive, as 



the elliptic functions must be evaluated at each integra- 
tion step, so the phase form of the equations is probably 
preferable if we are interested only in the exact solution 
to x{t). However, the averaged functions required for the 
adiabatic approximation to the solution are most easily 
derived from the analytic form of the equations, so this 
approach is better when we are interested in deriving an 
approximate solution to the equations. 



III. OSCULATING ELEMENTS FOR ORBITS IN 
THE KERR METRIC 

A. Summary of Notation 

In Boyer-Lindquist coordinates (<, r, 0, 0), the Kerr 
metric is 



2Mr\ , , Aasin^eMr , , 
1 ^ I dtdcl) 



2 „• 2 , 



+ (tu* - Aa'^ sin' 

S J 2 
+ irdr 
A 



Here 



E = + cos^ I 



A = r2 



a-" - 2Mr 



(3.1) 

(3.2a) 
(3.2b) 

(3.2c) 



and il/, a are the black hole mass and spin parameter. 
Throughout the rest of this paper we use units in which 
M = 1. for simplicity. 

We will make use of the Kinnersley null tetrad I, n, m, 
m* , which is given by 



I = + dr + , 



72 
2e' 



and 



1 



\/2{r + ia cos 1 _ 
The corresponding one-forms are 
1 = —dt + a siir^ 



Or H Oj, 

2S] ^^ 22 ^ 



ia sin 0dt + de 



(3.3) 
(3.4) 



smc 



A , aAsm^e 

n = dt H ( 

2E 2S 



-dr. 



dr 

2 



(3.5) 
(3.6) 

(3.7) 



and 



\/2(r + ia cos ( 



[—ia sin 6dt + T,d9 + im^ sin 



(3.8) 
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The basis vectors obey the orthonormahty relations l-n = 
— 1 and m-m* = 1, while all other inner products vanish. 
The metric can be written in terms of the basis one-forms 
as 

gai3 = -2/(an/3) + 2m(„m^) . (3.9) 
We define the conserved energy per unit rest mass /x: 

d 



dt 



(3.10) 



the conserved z-component of angular momentum di- 
vided by ^jlM: 



Lz=u - — , 



(3.11) 



and Carter constant divided by jjp'A'P: 

Q = ul- c? cos^ BE^ + cot^ BL\ + cos^ B . (3.12) 

(From now on we will for simplicity call these dimen- 
sionless quantities "energy," "angular momentum," and 
"Carter constant.") The geodesic equations can then be 
written in the form [iTj 



{Eir" + c?) - aL,^ - A [r^ + (L, - aEf + Q] 
Vr{r) , (3.13) 
Q - cot^ BLl ~ cos^ e{l ~ E^) 
VeiB) , (3.14) 



d<l) 
dX 



= csc^ + aE 



-1 - 



dt 
dX 



-TT = E 



V4r,B) , 

(r2 +a2)2 



- sin^ I 



aL-, 1 - 



(3.15) 

"2 + a^' 

A 
(3.16) 



Here A is the Mino time parameter [18[ , related to proper 
time T by 



d\ = —dr 

2-i 



(3.17) 



and we use these equations to define the potentials Vr(r), 
Ve{B), V^{r,B) and Vt{r,9). Sometimes it will be conve- 
nient to use, instead of the Carter constant Q, the quan- 
tity 



K^Q + {L,~ aEf 



(3.18) 



For the rest of this section we specialize to bound 
geodesies, which are periodic in r and B. 



B. Change of Variables 



Eqs. (|3.13p - p.l6p form a complete set of equations 
that can be solved to obtain the geodesic motion. How- 
ever, it is difficult in practice to use the variables r and 
B due to sign flips that occur in, for example, 

dr 



d\ 



at turning points. Therefore we follow Drasco and 
Hughes in switching to an alternative set of variables 



1. Angular motion 

We introduce the notation z = cos^ 0, and note that 
the effective potential can be written as 



Ve{z) 



1 



1-2 



[Q{l-z)-zLl~pz{l-z)\ , (3.19) 



where (3 = a^{l ~ E^). We note that this /3 is different 
from the variable appearing in the forced oscillator in 
Sec. HH All subsequent references to /3 will assume this 
new definition. We define Z- and 2;+ with z_ < 2+ to be 
the two roots, so that 



Veiz) 



1 



-/3{z- - z){z+ - z) 



(3.20) 



These roots z_ and z^ arc functions of E, and Q, and 
are positive with < z_ < 1 and z+ > 1. The motion 
takes place in the region < 2; < z_ . 

Wc replace the angular variable B, which oscillates, 
with another angular variable tpg, which increases mono- 
tonically with time. The definition is given by 



cos B 



z- cos ipg 



(3.21) 



Note that, for general forced motion z_, will change with 
time, along with B and ipe- 



2. Radial motion 

We define ri , r2 , and r4 to be the roots of the radial 
potential Vr{r): 

Vrir) = (1 - E^)in - r)(r - r2)(r - r3)(r - r^) . (3.22) 

Here the roots are ordered as < r4 < ra < 7*2 < ri, and 
the motion takes place in the region r2 < r < ri. The 
roots are functions of E, and Q. 

We replace the radial variable r, which oscillates, with 
an angular variable ipri which increases monotonically 
with time. The definition is given by 



P 



1 + e cos Ipr ' 



(3.23) 
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where the semilatus rectum p and eccentricity e are de- 
fined by 



ri = 



P 



1 - e 



r2 = 



P 



1 + e 



(3.24) 



where 



and 



n = L^-aE sin^ 9 , 



(3.31) 



C. Forced motion using tetrad components of 
acceleration and using convenient phase and energy 
coordinates on phase space 



We now turn to the forced geodesic equation 



dr^ 



(3.25) 



where a" is the external four acceleration. In this subsec- 
tion we derive our first formulation for integrating this 
equation, which parametrizes the acceleration in terms of 
its components on the Kinnersley null tetrad, and which 
parametrizes the motion in terms of a convenient set of 
coordinates on phase space that includes E, Lz and Q. 
The formulation is analogous to that presented in Sec. 
Ill Al for the nonlinear oscillator model. In particular, the 
phase variables used here are not conserved for geodesic 
motion. Our second formulation will be derived in the 
next subsection. 

Eqs. p.l3p - p.l6p are still valid for the forced geodesic 
equation. However, they must now be supplemented with 
evolution equations for E, Lz and Q (or K). We decom- 
pose the four acceleration on the Kinnersley tetrad as 



a ~ —anl — aifi + a*^rri -\- Gmm* 



(3.26) 



These four components are not all independent, since the 
acceleration must be orthogonal to the four velocity. We 
define 



Ra = 



and we take the three independent components to be a„, 
Ra and la- In Sec. lIII El we will show how the tetrad com- 
ponents of the acceleration, (a„, a;, a*^, a„i), relate to the 
acceleration components in Boyer-Lindquist coordinates, 
(at, ar,ag, a^). 

We similarly decompose the four velocity in terms of 
the Kinnersley tetrad as 



u = —Unl — uifi + it*„m -|- Umin* 
and we define 



(3.28) 



1 i 

Ru = —^('^"i + U*^) , /„ = -^{Ura - <„) , (3.29) 



The components are given by the expressions 

F 

A ' 
F A 



Ul — Ur 



21] 21] 



-Ur 



r oH cos 9 

Ru = ^U0 + . , 

2j Zj sm 9 

a cos 9 rTi. 
2^ Zj sm 9 



(3.30a) 
(3.30b) 
(3.30c) 
(3.30d) 



F = zu^E - aL, 



(3.32) 



The orthonormality condition u - a ~ allows us to solve 
for af. 

ai = - —an + — {RaRu + lalu) ■ (3.33) 

We also define the following three combinations of accel- 
eration components: 

Ai = rRa+ala cos 9, (3.34a) 
All = r la — aRa cos 9, (3.34b) 

Am = RuRa + lula ■ (3.34c) 

We can now write down the evolution equations for the 
energy, angular momentum and Carter constant. These 
are (see Appendix 



dE Urttn^ ^Aiii 



dX 



2Un 



a sin 9 All , (3.35) 



dLz a sin^ 9uran A asin^ 9AAin n 

—— = w sm9Aii 

d\ Un 2Un 



and 



dX 



(3.36) 
(3.37) 



1. Equations of motion in terms of phase variables 



We next replace the equations of motion p.l3p and 
p.l4p for r and 9 with new equations of motion for tpg and 
ipr, which are derived in Appendix The new equation 
for ipg is 



— = y^(I7^ 



1 



(1 — z^)Y,Ai cos^(Jg 



/3y^z^{z+ — Z-) sin ( 
cosipe sini/jgHaA^Aiii — 2wra„) 

2(z+ - Z-)(iUn 
cos ■061 sin iJeQAii 



P{z+-z-) 
where z = z- cos^ ijjg and 



(3.38) 



sin 9 



a^{l- Z-)sin9E 



(3.39) 
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The new equation for ipr is 



where 



dX 



CAiii sin tpr 



2(1 + e cos l/'r)Wn 2(1 + e COS tpryUn 



+ 



a£ sui9 simprAii 
1 + e cos ipr 

(l-e)2(l-cosV'r) 



Van 



u„(l + e COS^/jr)^ 



+ (1 + e)2(l + cosV'r) 



K2 



(3.40) 



7'=px/J/(l-e^) , 



(1 + e COS^r) 



il-E')^-^+[a'il-E') + Ll + Q] 



(1 + ecosV'r)^ , 



c = 

V 
Qi 



Qi(l-e) 22(1 + e) 



Kl K2 



5 = 



(1 - ey{l - cos^Ar)— + (1 + e)^(l + COSl/'r) 

Kl K2 

~2aL^rri - a'*£'(r + ri) + a^L^{r + ri) - a^S (r^ + r^ri + + rri(-2 + ri)) 

Erri [rri{r + ri) - 2 (r^ + rri + rl)) - [2a^E - 2Erri + aL^{-2 + r + ri)) cos^ 6 

-2aL,rr2 - a^E{r + + a^L,{r + r2) - a^E (r^ + r^ra + r| + rr2(-2 + 7-2)) 

i;rr2 (rr2(r + rj) - 2 (r^ + rrj + rl)) - (2a2£' - 2£;rr2 + aL^(-2 + r + ra)) cos^ 

i^i(l-e)(r + ri) F2(l + e)(r + ra) 



K2 



(3.41) 

(3.42) 
(3.43) 
(3.44) 

(3.45) 

(3.46) 
(3.47) 



r 



Here k = (r) , and subscripts 1 or 2 mean that a quan- 
tity is evaluated at r = ri or r = r2 (except for Qi and 
Q2). 



D. Forced motion using Boyer-Lindquist 
coordinate components of acceleration and phase 
variables that are conserved for geodesic motion 



In principle, we must evolve eight parameters, which 
are the four constants of motion and the four initial phase 
angles. However, one of these equations is eliminated by 
using the orthogonality condition x"aa = 0, where a dot 
denotes differentiation with respect to proper time, r, 
and a" is the acceleration. This condition is discussed 
in jlj and comes from the definition of proper time. 



In this subsection we derive our second formulation for 
integrating the forced geodesic equation, which is analo- 
gous to that presented in Sec. Ill Bl above for the nonhn- 
ear oscillator model. This formulation parametrizes the 
acceleration in terms of its Boyer-Lindquist coordinate 
components. It parametrizes the motion in terms of two 
phases "00 and xo defined below, which are conserved for 
geodesic motion, and three other parameters equivalent 
to E, Lz and Q, namely, the orbital eccentricity e, semi- 
latus rectum p and angle of inclination t (defined in Ref. 
[1^). This formulation is a generalization of the treat- 
ment of the Schwarzschild problem described by Pound 
and Poisson in [l|. 



In this section we will write the phase angles in the 
form "iAr = "0 ~ '00 J V'e = X ~ Xo a-nd derive explicit equa- 
tions for the time evolution of the initial-phase constants 
"00 and Xo- The other parts of the phases, ip and Xi are 
evolved using the standard geodesic expressions. While 
in practice we will need "0^ and "09 to evolve the orbit, 
we decompose the equations this way to facilitate com- 
parison to [l| and to make it easier to identify the con- 
servative contributions from the perturbing force, which 
are essentially ('0o)i (xo)- 
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1. Contravariant formulation 



The Gaussian perturbation equations (|1.3p carry over 
to the relativistic case and gives Eqs. (27)-(32) in [1|. In 
the Kerr case, we have two additional equations as the 
9 motion is no longer trivial. In [l|, the equations were 
integrated with respect to the anomaly. In the Kerr case, 
as we have two anomalies, it will be more convenient to 
integrate the equations with respect to the coordinate 
time, t. The equations of motion that are independent 
of the force terms are 

dr , dr I dr , dr , dr , 



do , do , do , do ,, 36 , 



(3.48) 
, 

(3.49) 



op oe Ob oyjQ oxo 



dt , dt . dt , dt . dt , . 

op Oe OL Oil-'o oxq 



(3.50) 
. 

(3.51) 



Here <I> and T denote the phase offsets for the evolution 
of (p and t. We can ignore these equations if we evolve t 
and (j) explicitly using the geodesic expressions evaluated 
along the instantaneous orbit, which amounts to evolving 
t — T and (/) — $ directly, as in the tetrad formulation. In 
the above, we use a dash to denote differentiation with 
respect to the parameter we use to define our orbit, which 
we take to be t. We will use a dot to denote differentiation 
with respect to the proper time r. The remaining four 
equations of motion are 



di I ^ di , ^ di I ^ di ^1 ^ di , 
dp^ de^ di dtpo ° dxo^^ 



dr I dr I dr , dr , dr , 
d-pP + d-e' + dl' + Wo^''' ^ 



dp^ de di dtpo " 9xo'^° 



dp^ de di dtpQ ° ^xo^^^ 



r / 

a T , 



The terms dr/dp etc. denote differentiation of the 
geodesic equations given earlier with respect to the var- 
ious orbital parameters. Following [H, we can use the 
orthogonality condition to get rid of one of these equa- 
tions, specifically Eq. p.52p . and we will directly inte- 
grate (j) and t which means we do not need to consider 
Eqs. (f330| - ([33T|) . 

We can rearrange Eqs. (|3.48p - (|3.49p to give 



dr f dr ^ dr ^ 
dr/dxpo \dp^ ^ de ^ di ' ' ^^^^^ 



de 



dO/dxo \dp 



do 



de 



(3.57) 



where we have made use of the fact that the equation 
for r, p.24p . is independent of xo and the equation for 6', 
p.2ip . is independent of -00 ■ The partial derivative dr/di 
also vanishes, but wc include this term explicitly for sim- 
plicity of notation in the following. Wc generalize [l| by 
writing 



Ca{x) 



dx 
da 



dr/da dx dO/da dx 



dr/dTpo di'o dO/dxo dxo 



(3.52) 
./ 

(3.53) 

r', 
(3.54) 

t' . 

(3.55) Substitution into Eqs. (|3.53l) - (|3.55p then gives 



(3.58) 
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p' = ^ {(Ce{e)C,{cf>) - CM^MV + {^dr)CM - ^^W^^^^^ , (3-59) 

e' = ^{{C,{e)Cp{c^) - C,{cf>)Cp{e))a- + {Cp{r)C,{cp) - C^^^^ , (3.60) 

L = ^ {{Cp{9)C,{cl,) - Cp{^)Ce{e))a'' + {C,{r)Cp{^) - Ce{c^^^^^ , (3.61) 

D = £p(r)(£e(e)A(0)-A(e)£e(0))-/:e(r)(/:p(0)A(0)- A(0)/:p(0)) + A(r)(/:j,(0)/:e(0)->Cp(0)/:e(0)) . 

(3.62) 



The correct evolution equations for the phase constants 
■00 and xo may be found by substituting the preceding 
equations into (|3.56p - p.57p . In the next section we will 
describe an alternative form of these equations which 
greatly simplifies the evolution of the constants of the 
motion. We include the above equations for completeness 
and to allow a direct comparison to the Schwarzschild re- 
sults described in [l|. 



through the constraint z^fa = 0. In the Kerr case, we 
do need to solve one of the radial or 9 equations, or some 
combination of them. Alternatively, using the definition 
of the Carter constant in terms of the Killing tensor, 
we can derive the evolution equation for Q straightfor- 
wardly. The time evolution of the related constant K 
defined in Eq. p.lSp . can be found from equation (|A9[) 
in appendix \K\ as K = K^^Uaap. The Killing tensor 
K"^ can be written in terms of I" and n" as 



2. Covariant formulation 



K 



(3.66) 



The preceding section presented the equations in a con- 
travariant formulation. We note that the equations for 
the evolution of the phase constants, p.56p - p.57|) . ap- 
pear to be singular at turning points where dr/dipo = 
or dd/dxo = 0. These are not real singularities, as the 
numerator also vanishes at the turning points, but it re- 
quires significant simplification to make this explicit. It 
is also possible to derive an alternative set of equations to 
p.52p - p.55p from a covariant formulation of the equa- 
tions. Pound and Poisson [l[ chose the contravariant for- 
mulation in the Schwarzschild case, since they found it 
easier to eliminate the singularities at turning points in 
that formulation. However, there are advantages to using 
the covariant formulation, since two of the covariant ve- 
locity components are then equal to conserved quantities, 
Ut — E, UAy — Lz- The osculation conditions become 



d^G fA 



, 



where 



9al3- 



dX 



(3.63) 



(3.64) 



in which /'^ denotes the orbital elements, including the 
phase constants. The first equation is the same as 
Eqs. (|3:48l) - (|33T|) which reduce to ^jSM^i an d (|33 7 |. Th e 
second equation is the equivalent of Eqs. p. 52^ - 1)3. 55p . 
but in this case two of the equations simplify significantly, 
namely, 



E^ft 



f<t> 



(3.65) 



In the Schwarzschild case, there is no equation for the 
6 motion and the radial equation follows from p.65|) 



from which we obtain 



K 



E—{w'^E~avj'^Lz)+Lz-^{a^Lz—a'CD'^E)—2Aurar , 

(3.67) 

where we have used E = —at-, and = a^. An alterna- 
tive expression for K"^ in terms of m" and m*" exists 
and is given in Appendix \X\ as Eq. (|A3[) . If we had used 
this definition we would have found an equivalent expres- 
sion for K that was a linear combination of E, Lz and 
ag. The two expressions are equivalent, since the orthog- 
onality relation between the perturbation force and four 
velocity always allows the elimination of one component 
of the force. 

These three equations provide an alternative way to 
evolve the constants of the motion, E, Lz and Q, but 
we must still evolve ipo and xo using p.56p - p.57p and 
therefore we still need to deal with the turning points. 

It is possible to derive an alternative form of these 
expressions that is manifestly finite at turning points by 
starting with the radial geodesic equation in the form 



T.^r^ = Vr{r,Lz,E,Q) 



We need to show that the term 



dr ■ dr ■ dr 
dE^^dLz^'^dQ^ 



(3.68) 



(3.69) 



that appears in an alternative version of Eq. p.56p . is 
proportional to r' . Differentiation of Eq. p.68p with re- 
spect to E yields 
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dr 
'dE 



dr 
dE 



2S2r — +21] 2r- 20^008 6* sin 



^d0_ 
'dE 



• 2 

r ~ 



dVr 

'dE 



dVr dr 
dr dE 



(3.70) 



Similar equations may be obtained by differentiating with respect to and Q. Multiplying the E equation by E etc. 
and adding the equations together, all terms on the left-hand side are proportional to r, while on the right-hand side 
we get the expression p.69p multiplied by dVr/dr plus the term 



dVr ■ , dVr ■ 

-dE^+dT^' 



A 

9g' 



= 2fE" 



1 dVr 

2S2"57 



(3.71) 



where the second equality follows from differentiation of Eq. (|3.68p with respect to time. The term in parentheses on 
the right-hand side is what we would obtain if we were on a geodesic, and therefore it necessarily equals in the 
evolving case. The final expression is 



dVr/dr 



■ dr 

^dE 



dr 
dlZ 



2^rrlE^ 



dr 

dL, 



do 



do 



A do 



2I]a^ cos 6* sin 6'f E—+L,^—— + Q— - S^a'' 



dE 



dL, 



dQ 



(3.72) 



in which V'gco denotes the geodesic expression for dtp / dr 
which we use to evolve ip- It is clear that this expres- 
sion is indeed finite at radial turning points, provided 
that the radial self-force is finite. In the Schwarzschild 
case, it may be easily verified that Eq. p.72p gives the 
same expression as Eq. p.56p when they are explicitly 
simplified. 

One important caveat is that although expres- 
sion (|3.72p is finite at radial turning points, it appears 
to diverge where dVr/dr = 0, and this condition will be 
satisfied once between each consecutive turning point. 
This is not a real divergence either, which is clear from 
the fact that the original form of the equations did not 
show such a divergence. Therefore, if we were to sub- 
stitute the various terms into the above expression we 
would find that the necessary cancellations would occur 
to eliminate these divergences. This simplification is a 
nontrivial calculation. However, an alternative approach 
that is easier to implement numerically is to use both 
Eqs. (|3.56p and p.72p without any attempt to simplify 
the expressions. By switching from one expression to the 
other near turning points we can avoid numerical round- 
off problems. This is the implementation that we use in 
practice and from which the results presented in Sec. IIVI 
were derived. We have verified in practice that both ex- 
pressions do yield the same results at points where nei- 
ther Vr nor dVr/dr vanish. 



action-angle formulation of the Kerr solution exists [19| . 
in which the equations take the form 

X = AxiE,L,,Q)Fx{^x-4'xo;E,L,,Q) , 

(3.73) 

dijjx 



dX 



= nx{E,L,,Q) 



(3.74) 



where X denotes (t, r, 9, </>) and A is "Mino time." The 
function Fx is periodic for r and 6, with a period of 2tt^ , 
and for t and it is the sum of a secular piece and an 
oscillatory term. The osculating element conditions give 



Fx 



Fx 



dA 



X 



F 



dE 

dAx 
dL, 

dAx 
dQ 



A 



dF 



X 



X 



A 



X 



A 



X- 



dE 

dFx 
dL, 
dFx 
dQ 



dE 
dA 

dA 
dQ 
dA 



d?A. 



xo 



dA 



(3.75) 



where the dash denotes differentiation of Fx with respect 
to the phase argument ipx — V'xo- As before, this ex- 
pression appears to be singular at turning points, where 
Fx ~ 0. However, we can obtain an alternative expres- 
sion by considering the potential 



3. Action-angle formulation 



dX 
dA 



— =Vx{X;E,L,,Q) . 



(3.76) 



The method described above for evolving the equa- 
tions of motion in the covariant formulation can be read- 
ily adapted to other problems and to other formula- 
tions of the Kerr geodesic solutions. In particular, an 



Adding the derivative of this expression with respect to 
E multiplied by dE/dX to the derivative with respect to 
Lz multiplied by dL^/dX and the derivative with respect 
to Q multiplied by dQ/dX gives 
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dVx 
dX 



dA 



X 



dE 



dFx 



dE 
dA 



F 



dA 



X 



X 



■A 



dFx\ dL 



X 



F 



dA 



X 



OF 



dA ■ V^'dQ 

dVx dE dVx dL 



Ax^]^ 



dQ 

dQ'J dA 

dVxdQ 



dX 



dA 



d fdX\ dE 
dE I dA j dA 



d 

dK. 



dE dA dL, dA dQ dA 
dX \ dL, _d_ /dX\ dQ' 
'dx)~dX^dQ \ dX ) dA 



The derivative of Eq. p.76p with respect to Mino time is 

dXd^X dVxdX dVxdE dVx dL, 



dA dA2 dX dA 
which thus allows us to replace Eq. p.75p with 



dVxdQ 

dE dA ~^ dL, dA ^ 9Q dA ' 



dVx . p, dV-A'o 



dX 



dA 



dX 
dX 



dV: 



X 



dX 



d^X 
'dA^ 



_d_ 

dE 



dX 
dA 



dE 
dA 



d 
dT, 



dX 
dA 



dL, 
dA 



_d_ 

dQ 



dX 
dA 



dQ 
dA 



(3.77) 



(3.78) 



(3.79) 



The term in square brackets vanishes for geodesies and is 
therefore proportional to the X component of the force 
when the orbit is perturbed. At turning points F^ and 
dX/dX are both zero and cancel, so we obtain a new 
form of the equation that is manifestly finite at turning 
points, albeit singular where dVx/dX = 0. As in the 
Boyer-Lindquist case, these two alternative formulations 
for the equations allow us to evolve the osculating ele- 
ment equations directly without worrying about singular 
behavior, just by switching between the two equivalent 
expressions in the vicinity of the turning points. 



E. Connection between Boyer-Lindquist and tetrad 
formulations 



The tetrad formulation of the osculation equations, de- 
scribed in Sec. IIII CI is written in terms of acceleration 
components, Ai etc., that are adapted to the Kinner- 
sley tetrad, while the Boyer-Lindquist coordinate for- 
mulation, described in Sec. IIIIDi is written in terms of 
the Boyer-Lindquist components of the acceleration. To 
identify the accelerations between the two approaches, 
we first write down the tetrad components of the acccl- 



^ The choice of periodicity is in a sense arbitrary, and different 
periodicities could be obtained by rescaling the angular variable 
ip. We specify a period of 2it for convenience. 



eration in terms of the Boyer-Lindquist components: 

(3.80) 
(3.81) 



tu^ A a 

a„ = — at ttr H ( 

21] 21] 2E 



ZD a 
ai = -^o.t + a. 



1 / i 

V2(r + ia cos 6) \ sin0 ^ ' 



(3.82) 



= —i^, 77 (-iasindat + ae 

(3.83) 

The acceleration functions Ra ~ {am + and la = 

i{a,n — o-m)/ introduced in Sec. IIII Cl havc components 



Ra 



a sin 9 cos 9 r a cot 9 



^ -at + —ae + ^ a^ , (3.84) 



ar sm U a cos u r 
la = — at - . a^ , (3.85) 

from which we obtain the tetrad acceleration components 
in terms of the Boyer-Lindquist components of the accel- 
eration 



All = —asiii9at 



2ra cos 9 1 



S sin6' ^ ' 

Am = —aa sin 9 at 

ug{r'^ — a^ cos^ 9)/ — 2ar cos 9 a 

H ae - ^—^H 

L, sm d 



(3.86) 
(3.87) 



in which 



a 



aE sin 9 ~ L, 
Esin6l 



(3.88) 



(3.89) 
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In Sec. IIVI below we will consider a toy problem as an il- 
lustration of the two methods. The force will be specified 
in Boyer-Lindquist coordinates, and the preceding ex- 
pressions can be used to obtain the corresponding tetrad 
components. 



F. Features and drawbacks of the two formulations 



By contrast, in the tetrad formulation used here, the 
magnitude of the acceleration is automatically finite. 
The independent components of the four acceleration are 
taken to be three of the four components on the Kinner- 
sley null tetrad, namely a„, a„i and a^, with the fourth 
component being determined by the orthogonality of the 
four acceleration and the four velocity. In terms of these 
three components, the square of the four acceleration is 



In this final subsection we discuss some of the advan- 
tages and disadvantages of our two formulations. 

First, as discussed in the Introduction, earlier work on 
methods of computing radiation reaction driven inspi- 
rals focused on the adiabatic limit [Tsl - fisl [20l . [2ll |. In 
this limit, it is sufficient to use orbit-averaged forces, 
or, equivalently, orbit-averaged proper time derivatives 
of the first integrals, E, and Q. These quantities can 
be computed as functions of i?, Lz and Q, both in post- 
Newtonian expansions and exactly using numerical black 
hole perturbation theory. In this paper our focus is on de- 
veloping methods that allow going beyond the adiabatic 
limit. For this purpose, orbit-averaged quantities are in- 
sufficient; one must use a prescription for the perturbing 
force that depends on the two nontrivial orbital phases. 
One could, in principle, continue to use the quantities E, 
Lz and Q to parametrize the force, if these quantities are 
taken to be functions of Lz and Q and of two addi- 
tional phases. This would be the most natural way to 
generalize the analyses of Refs. [isl - fTsi [20l[2l|. 

However, such a parametrization turns out to have a 
significant disadvantage compared to the parametriza- 
tions used in this paper, when one is attempting to com- 
pute approximate inspirals. Specifically, there are con- 
straints that the fluxes must satisfy at radial and polar 
turning points, in order to ensure that the four acceler- 
ation be flnite. Approximate versions of the fluxes may 
violate the constraints and lead to cusps in the motion at 
the turning points. (This will be true, in particular, for 
orbit-averaged fluxes.) The existence of these constraints 
can be seen from the expression for the square of the four 
acceleration in terms oi E, L, and K, which is 



1 



FF 



1 



5- 



gg 



(3.90) 



Here F = w^E ~ aLz, F = w^E - atz, g ^ asinOE - 
cscOLz, &ivdg = asinOE—cscOLz. It can be seen that at 
radial turning points where = 0, the fluxes must sat- 
isfy the constraint K = 2FF/A, while at polar turning 
points the constraint is _fir = 2gg. ^ 



^ A similar phenomenon occurred in the nonlinear oscillator model 
of Sec.|lT] where the time derivative of the energy was constrained 
to vanish at turning points. 



1 



(3.91) 

which is clearly always flnite.^ 

Similarly, in our Boyer-Lindquist formulation, the ac- 
celeration is again always finite, except in some special 
cases in the ergosphere. The independent components of 
the acceleration are taken to be the spatial, contravariant 
components a* = (a'', a^, a"^), with a* being determined 
by orthonormality. The square of the four acceleration is 
then 



9ij 



UiUj . 

gtt — o- I a a-' 



(3.92) 



which is always finite except in the ergosphere where it 
is possible for ut to vanish. 

We now turn to a discussion of a second issue, which 
is a significant advantage of the Boyer-Lindquist formu- 
lation over the tetrad formulation. This advantage is 
its simple behavior under the discrete symmetries of the 
Kerr spacetime. Specifically, note that any four acceler- 
ation a = a(a;",M^) can be uniquely decomposed as the 
sum of a dissipative piece and a conservative piece. For 
the dissipative piece, the components and are odd 
under u'' ~u^, — >■ ~u^, while the components a* 
and a'^ are even. For the conservative piece, the compo- 
nents a*" and are even, while the components a* and 
a'^ are odd [l2|. It follows that, in the Boyer-Lindquist 
formulation, wherein one specifies the components a*", 
and a*^ of the four acceleration, it is straightforward 
to independently specify the dissipative and conservative 
pieces. 

By contrast, in the tetrad formulation presented here, 
the independent variables are taken to be a„, am and a^, 
and the decomposition into conservative and dissipative 
pieces in terms of these variables is somewhat involved. 
In particular, if one is attempting to find useful approx- 
imations to the conservative self-force, for example, by 



^ A related issue is that the time derivative of the orbital eccen- 
tricity e can diverge oc 1/e as e — >■ 0, for forces parameterized in 
terms of E, Lz and Q, unless the fluxes obey certain constraints 
at e = 0. This issue is discussed in detail in Ref. [13 . Again, 
this divergence is automatically excluded if one parameterizes the 
force in terms of its tetrad components; The eccentricity can be 
written as a smooth function e = e{x°' ,p^) on phase space. Tak- 
ing a proper time derivative gives de/dr = ma°'de/dp°', which 
is finite for finite accelerations. 
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naively using conservative post-Newtonian approxima- 
tions to the quantities a„, and a'^, the errors in the 
approximation will generically lead to a self-force with 
both conservative and dissipative pieces. This can be a 
problem since in the adiabatic limit the effect of the dis- 
sipative self-force on the motion is boosted relative to the 
conservative self-force. 

There are alternative parametrizations of the self-force 
that combine the advantages of our two formulations, for 
example, 

a" = a^ef + Je^ + a^e^'^^gu'^ejel + {a^Uf + Jug)u", 

(3.93) 

where eV and Cg arc unit vectors in the directions of dr 
and de. Here the dissipative and conservative pieces of 
the quantities a^, and a± have simple transformation 
properties under discrete symmetries, and moreover the 
magnitude of the four acceleration is 



{ar [1 



+ ia' 



e\2 



1 



2 I 2 



(3.94) 

which is always finite. Useful approximation schemes can 
be obtained by (i) formulating approximations in terms 
of the three variables a^, and a±; (ii) using the exact, 
Kerr relations to compute a„, am and in terms of 
a^, and a±; and (iii) using the resulting expressions in 
the tetrad formulation equations of motion (|3.15l) , p.l6p , 
(1?:^ - ^Ml, (IX^ and (IXiOl) . See Ref. for an 
application of this approach. 



IV. EXAMPLE OF PERTURBED KERR 
GEODESICS: "GAS-DRAG" 

As an example problem, we will suppose that the small 
mass experiences a drag force proportional to velocity, 
which could represent the behavior of an EMRI occur- 
ring in the presence of gas. Here we derive the four ac- 
celeration for such a force. 

In a given frame of reference, the relativistic analog of 
this simple drag force will have a term proportional to the 
spatial part of the velocity, plus a term proportional to 
the frame velocity constructed so that the force remains 
orthogonal to the total velocity. Let wzamo be the veloc- 
ity of zcro-angular-momentum observers (ZAMOs), and 
let u be the velocity of the small mass. In the frame of a 
ZAMO, the spatial part of the velocity of the small mass 
is 

u± =u + rwzAMO , (4.1) 

where T = u ■ wzamo- The drag force then has the form 

a = -7Uj^ + KUzAMO -ju + {k~ jT)uzamo , (4.2) 

where 7 is the linear drag coefficient. Enforcing the con- 
dition a ■ il = then determines the value of k 



7(r2 - 1) 



Inserting this into the formula for the acceleration due to 
drag (|4.2|) gives 



a = —7 u 



U ■ ItZAMO 



(4.4) 



Writing this explicitly in terms of Boyer-Lindquist coor- 
dinates, we have 



(4.5) 



a = —7 u 



where u° denotes the four velocity of the inspiraling ob- 
ject and 



'(r2 -|-a2)2 - Aa2 sin^ 6* 
2ar 

SA((r2 + a2)2 ~ Aa2 sin^ ( 



2r \ J 2a sin Or 



'■z ' 



(4.6) 
(4.7) 

(4.8) 
(4.9) 



in which A 



S ~ ^a} cos^ d as before. 



(4.3) 



As a test case, we constructed an inspiral into a central 
black hole with spin a = 0.9 under the influence of this 
gas-drag force with 7 = 10^''. Wc took the initial orbital 
parameters to be vjM = 7, e = 0.5, i = 7r/6, = t = 0, 
■0^ = 1 and "06 = 2. The inspiral trajectory was con- 
structed using both the Boyer-Lindquist and the tetrad 
formulations. The evolutions were found to be identical, 
as we would hope, and this gives us confidence that our 
results are correct. In the following discussion, we will 
not distinguish between the results obtained using the 
different formulations as they differed only at the level of 
numerical noise. 

In Fig. [2] wc show the evolution of the orbit under 
the infiuence of the gas-drag force and initial condi- 
tions given above. The six panels show the three Boyer- 
Lindquist coordinates, (r, 0, 0), and the three constants 
that describe the orbital shape, (p, e, t), as functions of 
the Boyer-Lindquist time t. We see that the influence of 
the drag force is to drive the inspiral of the object, but 
also to increase the eccentricity of the orbit and decrease 
the orbital inclination, i.e., to make the orbit more pro- 
grade. The trajectories of the geodesic constants of mo- 
tion, (p, e, t), show oscillations on the orbital time scale, 
superimposed on a monotonic secular evolution on the 
radiation reaction time scale. The secular part is the 
analog of the averaged evolution, G'^^-', described for the 
perturbed nonlinear oscillator in Sec. [Hi The panels in 
Fig. [5] also show the solution to the adiabatic equations 
of motion for this problem computed using Eqs. (|2.8p - 
(|2.15p . We see that the adiabatic solution is a good ap- 
proximation to the average evolution along the inspiral, 
as expected, and, as we saw for the toy problem in Sec.llll 
it provides a closer fit to the change in the constants of 
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FIG. 2: Evolution of the orbit under the influence of the "gas-drag" force. The panels show, as functions of Boyer-Lindquist 
time t, the particle coordinates r (top left), 6 (top right) and (p (middle left) and the orbital constants p (middle right), e 
(bottom left) and l (bottom right). In each panel the solid curve was computed using the exact evolution, while the dashed 
curve was computed using the adiabatic evolution. The plots showing (j), p, and l have insets which show close-up views of the 
same data so that the difference between the adiabatic and exact results can be seen. 
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the motion, in this case p, e and than for the phase. 
In this case, although the adiabatic solution remains in 
phase on average over the whole of the evolution seen 
in Fig. [5J within each cycle the adiabatic solution goes 
in and out of phase with the true evolution. This arises 
because the forcing term in this case is relatively large, 
and so the orbit changes significantly between periapsc 
and apoapse as a result of the forcing term. The or- 
bit evolved using the instantaneous force therefore looks 
quite different from the orbit evolved continuously by the 
orbital-averaged force. Note that the orbits do come back 
into phase after each complete cycle, as expected. 

This example illustrates the application of the oscu- 
lating elements formalism to the computation of inspi- 
ral evolutions in the Kerr spacetime, and it serves to 
demonstrate that the two alternative formulations do in- 
deed yield the same results. However, even though the 
prescription for the drag force was rather simple. Fig. [5] 
also illustrates some qualitative effects of the drag force 
that could be used to infer the presence of such a drag 
from observations. For an orbit evolving under the in- 
fluence of gravitational radiation reaction only, the ec- 
centricity tends to decrease, except toward the end of 
the inspiral just prior to plunge [2^ [l^, while the in- 
clination tends to increase, i.e., the orbit becomes more 
retrograde [13, HH, HJ] . We see here that the effect of the 
drag force is qualitatively different, as it drives increas- 
ing eccentricity and decreasing inclination. If observed, 
this would provide a robust observational signature for 
an orbit that was evolving under the influence of drag. 
A decrease in orbital inclination due to hydrodynamic 
drag was also seen in [25| , in which a more sophisticated 
model for the drag force was employed. It is gratifying 
that this simple model produces this expected feature 
qualitatively. The same paper [25| found that the eccen- 
tricity would increase in parts of the parameter space and 
decrease in other parts. For Newtonian orbits, the oscu- 
lating element equations predict that the eccentricity will 
remain constant under the action of a simple drag force of 
this type (see Appendix [C|) . Increasing eccentricity has, 
however, been seen in Newtonian simulations of binaries 
embedded in a realistic disc f26l-l28j. In Appendix ID] we 
show that the increase in eccentricity is an expected post- 
Newtonian effect and give an explanation in the context 
of a Schwarzschild black hole. It is clear from that discus- 
sion that this increase in eccentricity is a generic feature 
of relativistic drag, and so this is an observational predic- 
tion. If observed, an increasing eccentricity or decreasing 
inclination would be a clear signature that the observed 
inspiral was not occurring in a vacuum Kerr background. 



V. DISCUSSION 

We have described two methods for integrating the 
equations of motion for bound, accelerated orbits in the 
Kerr spacetime, which are based on identifying the or- 
bit with a geodesic at each point. The first method 



parametrizes the position and velocity of the orbit in 
terms of the conserved quantities (energy, axial angu- 
lar momentum and Carter constant) in addition to three 
angular variables which increase monotonically and cor- 
respond to relativistic generalizations of the anomalies of 
Kcplerian motion. The second method is the traditional 
"osculating element" technique which parametrizes the 
position and velocity of the orbit in terms of the geodesic 
with the same position and velocity. Practically, the sec- 
ond method differs from the first only in the treatment 
of the three phase variables, which are split up into a 
geodesic piece and a "phase offset" piece that is constant 
for geodesies. 

To illustrate the methods, we first analyzed, as a sim- 
pler model, a forced anharmonic oscillator. This was 
written in terms of a set of phase space coordinates. The 
forced equations of motion contained an apparent diver- 
gence at the turning points, but it was possible to refor- 
mulate the equations to eliminate the problematic terms 
and thus obtain equations of motion in a form without 
divergences. We discussed the adiabatic prescription for 
computing the leading order motion, which corresponds 
to a gradual evolution of the oscillator's amplitude and 
fundamental frequency driven by the phase space aver- 
aged forcing function for the amplitude. We presented 
an alternative analysis of this toy problem analogous to 
the osculating orbit method in terms of the analytic solu- 
tion to the un-forced motion. By numerically integrating 
the equations we verified that both parametrizations gave 
the same results and compared these to the adiabatic ap- 
proximation to the solution. 

Next, we showed that the equation of forced motion 
in the Kerr spacetime could be reformulated in a simi- 
lar fashion. For the first method, it was advantageous 
to parametrize the force in terms of its components on 
the Kinncrsley tetrad instead of using the instantaneous 
time derivatives of the conserved quantities. We derived 
a formulation of the equations of motion in terms of phase 
variables that was manifestly divergence-free at the turn- 
ing points. We then generalized the second method, of 
osculating orbits, to generic orbits in the Kerr spacetime 
and showed how we could write down a divergent-free 
form of equations of this type without explicit simplifi- 
cation. 

As an application of our results, we considered the case 
of a simple force that could represent a gas drag. Numer- 
ical integrations of the equations of motion for a choice of 
parameters verified that the two methods of parametriz- 
ing the motion gave the same results. We identified a key 
observational signature of the presence of a drag force, 
namely, a decrease in the orbital inclination and an in- 
crease in eccentricity, which is opposite to the increase 
in inclination and decrease in eccentricity characteristic 
of the gravitational radiation reaction forces during the 
early stage of an inspiral. 

The first of our two methods has been applied to the 
study of transient resonances that occur in the radiation- 
rcaction-drivcn inspirals of point particles into spinning 
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black holes, using approximate post-Newtonian expres- 
sions for the self- force [22j. Other applications of this 
work will include the construction of accurate trajecto- 
ries for orbits evolving under the action of the self-force, 
once self-force data for generic orbits are available. This 
will be essential for the construction of accurate gravi- 
tational waveforms for EMRIs, which will be needed for 
LISA data analysis. The formalism can also be used to 
estimate the magnitude of any secular changes in the or- 
bital parameters that arise from the action of external 
perturbing forces. These could arise from gravitational 
perturbations from distant objects, such as stars or a 
second massive black hole, or from the presence of other 
material in the spacetime, such as the gas-drag which we 
considered in a simple way here. It will be very important 
to have a quantitative understanding of the importance 
of all these effects if intcrmediate-mass-ratio inspirals or 
EMRIs are to be used to carry out high-precision map- 
ping of the spacetime around Kerr black holes and for 
tests of general relativity. Finally, the results described 
here will be useful to augment existing kludge models 
for inspiral waveforms. In particular, these methods will 
allow us to extract the secular part of the evolution of 
both the orbital constants of the motion and the phase 
constants, from self-force calculations. It is straightfor- 
ward to include secular changes to the orbital parame- 
ters in the kludge framework [l3| . and by doing this it 
should be possible to ensure that the kludge waveform 
stays in phase with the true waveform for long stretches 
of the inspiral. It will be important to have accurate but 
cheap-to-calculate waveform models available when the 
data from gravitational wave detectors arc analyzed, as 
this data analysis will rely heavily on matched filtering 
using template waveforms. 
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Appendix A: Derivation of tetrad equations of 
motion in terms of radial and polar angular variables 



In this appendix we derive the forms Eqs. p.35p - 
p.37|) . p.38p and (|3.40p of the equations of motion for 
forced motion in Kerr, in the tetrad formulation, and 
using the angular variables V'r and ipe instead of r and 6. 



1. Evolution equations for first integrals E, Lz, K 



The evolution equations p.35p - p.37p for the con- 
served quantities E, Lz and K are obtained as follows. 
We start from the standard expressions for the first inte- 
grals in terms of the Killing vectors and Killing tensors 

r = -sr , (Ai) 

^ S^}, (A2) 
= 2Em("m*'') - cos^ 6*5"'^ , (A3) 

and take proper time derivatives. Using the tetrad de- 
composition (|3.26p of the acceleration together with the 
expressions (|3.6p - p.Sp for the basis covectors gives 
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Here, we have used the fact that m given in Eq. p.8p can 
be written as 

1 



m = ^ [—{ir + a cos 9) sin 9dt + (r — ia cos 0)T,d9 

(A5) 
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Noting that «,„ ~ {Ra~ila)/ \/2, and eliminating ai with 
the aid of Eq. p.33p transforms Eq. (|X4)) to the form 



dE 
rf7 



a„ / A 

Un - —Ui 

Un \ 2L 



asmt 



A 



2Sw„ 



{RaRu + lalu) 



{ria — a COS 9 Ra) 



(A6) 



Using the expressions p.30p for the tetrad components 
Un and ui of the four velocity and converting from t 
derivatives to A derivatives using Eq. p.l7p then leads to 
the final form given in Eq. p.35p . 

Similarly, we obtain Eq. p.36p by starting from 
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and eliminating a/ using Eq. p.33|) to obtain 
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The form quoted in Eq. p.36p is then obtained from this 
using Eqs. (|3.30p as before. 

The evohition of the Carter constant is obtained very 
simply from the expression for the Killing tensor to be 

^ = 2K''^u^afi = 2S (RuRa + luh) , (A9) 

(XT 

where we have used the orthogonality relation g'^^Uaap = 
and written the combination {u^am + Umo!^) in terms 
of Ru and Ra- Combining this with the definition (|3.34c|) 
oiAiii yields Eq. (|337l) . 



Polar motion 



To obtain the equation of motion p.38p for if^g , we start 
by differentiating its definition (|3.2ip with respect to A: 
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The equation of motion p.l4p for 6 can be rewritten in 
the form 



d9\^ , ^ (z+ - cos^ ipe) 

— = Pz^sm iPe--. 

dX J (1 — z_ cos'^ ■i/'e) 



(All) 



We now take the square root of this equation. From the 
definition p.2ip of "tpg and noting that ijjg monotonically 
increases, we see that [dO/dX) > Q ioi Q < ipg < ii and 
[dO/dX) < for TT < -0e < 27r, so we must choose the 
positive square root on the right-hand side. Combining 
Eq. (|A10p and the square root of Eq. (|Alip now leads to 
an equation of motion for tpg in the form 



— = V/?(^+-z.cos2^,) + ^ — . (A12) 



Next, we can obtain an expression for dz- / dX in terms 
of dPi/dr, where Pi = {E^L^.K) are the constants of 
motion, as follows. Using the chain rule gives dz^/dX = 
{dz-/dPi){dPi/dX). Differentiating Vg in the form given 
in Eq. p.20p with respect to Pi at fixed z and evaluating 
the result at z_ relates dz^/dPi to {dVg/dPi)z_ , which 
can be computed from Eq. (|3.19p . This yields 
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Now switching from the Carter constant Q to K = Q 
[Lz — aE)'^ and using that dX — dr/T, we obtain 
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The expressions for the evolution of Pi in Eqs. (|3.35p - 
p.37|) can now be used to obtain the explicit dependence 
on ipg of dz-/dX given by Eq. (jA14p by direct substitu- 
tion. This gives 



dz- 2A'H^Ura„ 

dX Uri 



— 2_)— ^ = — — ^ az_ sin^ -06 + 



2(1 - 2:_)(r^ + a^z_ cos^ Ve) - az^y.- 



A sin ^g 



+271- \/\ — z- cos^ (^^ + a^z-) (rla — acosdRa) 
where H- = H{z-) = — a(l — z^)E. This can be written as 
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Substituting the expressions for the four velocity components Eqs. p.30p and the definition p. 311) of Ti into the 
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coefficients of Ra and la inside the square brackets and expanding them out gives 
dz A'H 

— — — az- sin^ 'ipg [Am — 2uran] + 2{1 — z-)'E,ue{acos9Ia + rRa) 

+ -^Z=^£^[vD^L,-a^E{l- z^){l- z. cos'' ije)]{rla- acos0Ra) . (A17) 

^1 — COS^ ^0 

Using that ug — (dd/dX), with (dd/dX) given by the positive square root of Eq. (|A11|) and inserting Eq. (|A17P into 
the equation of motion for i/jg of Eq. ()A12|) leads to the fina l resuh quoted in Eq. p.38|). 



Radial motion 



We now give a derivation of the radial equation of mo- 
tion p.40p which is similar to the above derivation of the 
equation p.38p of polar motion. From the definitions 
p.l3|) and (|3.22p of the radial potential we have 
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where F was defined in Eq. p.32p . We parametrize the 
roots of the right-hand side by Eq. p.24p for the turning 
points ri and r2 of the bound motion, and by 
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for the other two roots. Substituting the definition p.23p 
of ipr into Eq. (|M8)) and using Eqs. ([XM]) and (|M9)) 

gives, after some algebra, 
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By differentiating the definition p.23p of ipr we obtain 
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We note that ipr is chosen to monotonically increase, 
which means dipr/dX > 0. We specialize to the con- 
vention that ■i/^r = at r = r2 and ipr = tt at r = ri, 
so that r increases tor < ipr < and decreases for 
TT < -0^ < 27r and we choose the positive square root 
i n Eq . ((MO)) . Substituting Eq. (|M0)) for dr/dX in Eq. 
(|A2ip shows that the geodesic term becomes 



dipr 
'dX 



[p{l ~ e) ~P3{1 + ecosV'r)] 



1/2 



n _ e2 , 

ycodcsic V / 

1 /2 

X + e) - p4(l-|- ecos-^r)] ' 

= V . (A22) 



Here one can check that Eq. (jA22[) is just a 
reparametrization of Eq. (j3.42|) by substituting the ra- 
dial potential in the form given in Eq. p.l3p in terms of 
P, into Eq. dMU, 

since V = (dr/dipr) ^sfVr-, expressed 

in terms of 'ip^- 

The nongeodesic terms in Eq. (jA21[) are obtained as 
follows. From Eq. (|3.24p for r\ and r-i it follows that 
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Substituting Eqs. (|A23p into Eq. (|A21|) gives 
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Next, expressions for the derivatives of the turning 
points ri and r2 can be computed in terms of dPi/dX 
by using that {dri^2/dX) ~ {dri^2/dPi)dPi/dX. Differen- 
tiating the radial potential with respect to Pi at fixed r 
and evaluating the result at ri and r2 gives 
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= -(1 - E')in - r2){r2 - r,){r2 - r4)|g . 

(A26) 



We note that one can see from Eqs. p.22p . (|A18|) and 
(|M5)) - (U26|) that the coefficients of dri,2/dP^ can be ex- 
pressed in terms of the r-derivative of Vr at fixed Pi eval- 
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uated at the turning points as 
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Here K(r) = V^{r), which can be computed from Eq. 
ens)) to be 



where the definition (|3.32p of F has been used. Using the 
derivatives of Eq. (jAlSp with respect to Pi then results 
in the following expressions for dri^2/dX: 
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K(r) = AEFr - 2rA - 2(r - M){r^ + K) , (A29) With this, Eq. (|MT|) becomes 



dTjjr 



1 



(1 - e)^{coa'ipr 



2ep sin -0^ 

-(1 + ef{cos7pr + 1) 



1) 


' 


1 ZU 


Kl \ 


2F2 


f 2dE 


K2 





^dE_ 
^ dX 

dL, 



dL 



Ai dK 
dX J Kl dX 
A2 dK' 
K2 dX 



(A31) 



The next step is to substitute the expressions p.35p - p. 371) for the derivatives of the first integrals into Eq. (jA31l) . 
After some algebra we obtain 
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where Si = mj-a'^ sin^ 9. Noting that Ur = A^^{dr/dX) Substitution of Eq. (|A33|) and Eq. ((3?30)) together with 

and using the definition (|A22p of P gives an explicit ex- further algebraic manipulations on Eq. (|A32p lead to 
pression for u^: 
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We can simplify the coefficients of Ra and la by expanding the term {1 
in Eq. p.2p to obtain an explicit factor of [r — ri): 
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where z — cos^ 6, as before. We similarly define Q2 by replacing 1 ^ 2 in Eq. (jA37|l . Substituting Eqs. (|A37|) 
as Eqs. (|A33[) and (jA35[) into Eq. (|A36p and using the definitions p.34p yields after simplifications 
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This can be further simplified to be 
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Appendix B: Adiabatic Limit 



In this appendix, we derive our method of obtain- 
ing the leading order, adiabatic solutions to the forced 
geodesic equations in Kerr. This method was used to 
obtain the numerical adiabatic solutions that are plotted 
and discussed in Sec. II VI above. The starting point is the 
specific form p.35p - (|3.40p of the forced geodesic equa- 
tions derived in Sec. IIIIC] above, which have the general 



form 

V'a = c^a(V'a,J) + e5i'^(^,J) + 0(e2), 

1 < a < AT, (Bla) 

Ja = eGi')(t/>,J) + e2Gf (V',J) + 0(e3), 

1 < A < M. (Bib) 

Here xjj = {ipi, . . . are a set of angular variables, 

and J = (Ji, . . . , Jm) are a set of quantities that are con- 
served for the unperturbed system. Dots denote deriva- 
tives with respect to A. The functions determine the 
frequencies of the unperturbed motion (geodesic motion 
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for the Kerr application) , and the functions ga , Gj^ and 

G^^' represent the external perturbations on the system*. 
These functions are all periodic in each phase variable 
with period 27r. In the special case when the frequen- 
cies Wa are independent of the phase variables ip, the 
variables ^pa and Jx are (generalized versions of) action- 
angle variables. This special case is actually fully general; 
one can always perform a redefinition of the phase vari- 
ables to achieve this. This case of action-angle variables 
was studied in detail in Ref. [l^l, where the form of the 
adiabatic and post-adiabatic solutions were derived. 

Here we will generalize the analysis of Ref. [l^l to the 
more general system of Eqs. (|Bip . since our system of 
Eqs. in Kerr is of this form. We start by 

describing the result for the adiabatic limit, and then we 
outline its derivation. The adiabatic solutions are given 
by the following set of steps: 

1. We define the averaging operation, for any function 
f{ip) of if), by 



Jo ■ • ■ Jo t^iv(V)jv,j) J^yi' ■ ■ • 'y^-* 
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The subscript J on the left-hand side is a reminder 
that the averaging operation depends on the value 
of J. 

2. We define the averaged frequencies and forcing 
functions 
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3. We solve a set of ordinary differential equations in 
the slow time parameter 



A = eA 
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for two sets of auxiliary functions and J^\{X). 

This set of ordinary differential equations is 
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Note that for this step, one does not need to specify 
a value of e. 



* Note that the notation u)a{'4>a,^) means that each u)a depends 
only on a single phase variable Va, and does not depend on 
the phase variables 'ijip with fi ^ a. The adiabatic limit of the 
more general system of equations with uja = tjJa{ip,3) would be 
considerably more complicated. 



4. We can then write down the adiabatic solutions: 

(B7a) 
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where the function ^a{Xi J) is defined implicitly by 
the equation 
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S„(x + 27r,J) = S„(x,J) + 27r 



(B8) 
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We now turn to the derivation of this result. We start 
by rewriting the differential Eqs. (jBip in terms of the new 
variables (x^, Jx), defined implicitly by the relation 



■0q(Xq, J) = ^^(Xa, J) 



(BIO) 



All of the functions appearing in the differential equa- 
tions are expressed as functions of the new phases Xa i 
they must be periodic functions of each Xa by virtue of 
the property (lB9|. Using the definitions (HH, (|B2]) and 
(IB3I) the result can be written in the form 



Jx 



-4J)+.^^^^5i^)(x,J) + 0(e^) 
l<a<N, 

eGi^' (X, J) + e'G^x^ {x, J) + 0{e^), 1<X<M 



(Blla) 
< M . 
(Bllb) 



This system of equations is now in a form to which the 
results of Ref. [l2| can be applied; the variables (xa, Jx) 
are generalized action-angle variables. The averaging op- 
eration defined in [l2j . a straightforward averaging with 
respect to the phases Xqi coincides with the definition 
(jB2p used here, because of the definition (jBlOl) . The re- 
sults of Ref. [l2| now imply that the leading order solu- 
tion for Jx is of the form given by Eqs. (|B7a| and (|B6b|. 
They also imply that the leading order solution for Xa 
is of the form Xq(A, e) = Xa(eA)/e, where Xa satisfies 
the differential equation (jB6ap . Combining this with the 
definition (|B10p now yields the result (|B7b[) . 



Appendix C: Perturbation of Keplerian Orbits 

Here we derive the osculating element equations for a 
Keplerian orbit experiencing a force in the plane of the 
orbit, f = — /.tr/r'^. In this case, we can take the orbital 
plane to be the x-y plane. The orbit is described by four 
parameters — the semimajor axis, a, the eccentricity, e, 
the argument of perihelion, w, and the time of pericenter 
passage, Tq. (The restriction to a plane gets rid of the 
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other two orbital constants.) The orbit is elliptical and 
described by 

r = ^ail-ecosE) , (CI) 

l + ecos(u-w) ^ ; ' V ; 

^3(1^^2)3 (1 + e C0S(W - W))' , (C2) 

in which u is the argument. It is usual to call v = u — uj 
the true anomaly and E defined by the first equation 
above is the eccentric anomaly. The time of pericenter 
passage is given implicitly by 



dv' 



Q (l + ecosw')^ Y a'^(l — e^)'^ 



{to -To) , (C3) 



where wq = w(io)- 

Under the action of a force in the orbital plane with 
radial component R' and tangential component S", the 
Gaussian perturbation equations predict the following 
evolution equations for the four orbital elements [lo| : 



'a(l-e2) 2a 
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e = J ° ^ ^ ^ ^ [sin vR' + (cos v + cos E)S'] , (C5) 
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If we consider the true anomaly, v, then since i; = u — w, 
V = u — ui. By the definition of the osculating elements, 
the value of ii is always given by the geodesic value, 
and so we see that the evolution of the true anomaly 
differs from integrating the instantaneous time-evolving 
geodesic equation by the w term. This can also be seen 
by differentiating the orbit equation and using that both 
r and f are consistent with the instantaneous geodesic to 
obtain 



y/fiajl - e^) 



e cosu 



e smv 
2ee \ / 1 + e cos v 
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The first term is the geodesic v, while the other terms 
arise as a result of the perturbation. Although this equa- 



tion looks singular at turning points, smv = 0, substi- 
tution of the expressions for a, e and the geodesic equa- 
tions gives the necessary calculations and the expression 
reduces to — cj, as it should. 

Appendix D: Drag force in Schwarzschild geometry 

In order to understand the effect that leads to an in- 
crease of eccentricity we can consider a Schwarzschild BH 
system, in which the same effect is seen, but which is eas- 
ier to analyze and to understand. The osculating element 
equation for the evolution of the eccentricity, Eq. p.60t . 
in the case of a nonrotating BH reduces to Eq. (37) in [l| 
and has the form 



de 
dv 



TZ{p, e, v)a^ + Tip, e, v)a'^ 



(Dl) 



We use a drag force to perturb the orbit which takes a 
very simple form a"^ ~ —ju^, a"^ = —ju'^. The velocities, 
in Schwarzschild coordinates, are 



p — 6 — 2e cosv 



V pip — 3 — e^) 
(1 + e cosw)^ 



(D2) 
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The equation for de/dv is integrable for this perturb- 
ing force if changes to e and p are ignored over the orbit 
and the result can expressed in terms of elliptic integrals. 
However, this is quite messy and we are primarily inter- 
ested in the leading order correction to the orbit. We 
make a weak field expansion [M/p « 1) of the terms 
entering this equation: 



ff (^o(e,^;) + y7^l(e,«) + 0(^fVp2)^ ^ ^4) 
r « pHl v) + ^Ti{e, v) + 0(MV/)) , (D5) 

here we do not go beyond the first correction to the Ke- 
plerian term. Similarly, we find for the velocities 



ur = ^l^ul,{e,v) (^1 + jul{e,v) + 0{M^ /p^) ) , 



n^^J ^4{e, v) (l + ^u?(e, v) + OiMyp') 



(D6) 



(D7) 



The explicit form of the terms in these expansions is 
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To 
Ti 



e smf , u[ 



e cos 7J H — e' 

2 2 



'3 1 

(1 + ecos^;)^ = ( ^ + 



7^l = 37eo(l - e^) 



(1 + ecosi;)^ 
(ecosw + 2) cosw + e 

(1 + ecosu)^ ' 
2e + 6 cos V + 2e cos^ v — e'^ cos'^ w — 5e^ cos v — cos^ v — 
(1 + e cosw)* 

I 



(D8) 
(D9) 
(DIO) 
(Dll) 
(D12) 



The leading order terms give us the Newtonian perturba- 
tion of the eccentricity (|C5[) with perturbing force com- 
ponents R' ~ ~jr, S' = —^rcf). Overah, the Newtonian 
term is 



-2-fp 



3/2. 



cos V 



00 



(1 + e cosf )^ 



(D13) 



This equation can be integrated over an orbit, keeping 
e,p on the right-hand side constant, to give 



Se(v) 



^2-fp 



3/2_ 



smw 



1 



e cos V 



(D14) 



It is clear that in the Newtonian case there is no secular 
change in the eccentricity. Note also that the individual 
components (radial and azimuthal) of the perturbation 
are not zero after integration over one orbit, but they 
are exactly equal and opposite in sign. We now consider 
the first relativistic corrections. First, we note that the 
perturbations TZi and uf are independent of v, and so we 
can reabsorb these into a redefinition of 7 — > 7' where 



7 = 



1 



P 



7 



and so the rescaled leading order term still averages to 
zero, as it is proportional to the Newtonian expression. 
There remain two perturbations, one that comes from the 
radial velocity perturbation, TZqu\, and one that comes 
from the relativistic correction to the orbit's response to 
the azimuthal perturbation, u^Ti. The velocity pertur- 
bation contributes 



01 



1/9 -2 

, p ' e sm V 
(1 + ecosu)^ 

-(-3 - 2ecosw + e^) - -{3 + e^) 



^ ^ (1 + e cos vY 



(D15) 



Note that this term is always positive and so it will lead 
to an increase in the eccentricity. This can be interpreted 



as an additional radial force which acts at each point of 
the orbit in the direction of motion slowing down the 
effective radial velocity in the force, which leads to the 
increase of eccentricity. 

The second part of the perturbation, Uq7i, contributes 



7'epi/2 



X [1 — + cos^ v[l 4- ecosw) — e(cosw + e)] 
X (1 -f ecosu)-2 . (D16) 

Note that the last term is proportional to the Newtonian 
term and therefore averages to zero so we can ignore this 
term. The remaining part is always positive and also 
drives an increase in eccentricity. This time the extra 
term can be interpreted as an additional azimuthal force 
which further boosts the effective azimuthal velocity in 
the force and once again leads to the increase in the ec- 
centricity. 

We note that both of these perturbations, and also 
the Keplerian term, are proportional to eccentricity, and 
they will not drive a circular orbit to become eccentric. 
In fact, the contribution from the relativistic correction 
to the velocity is equal to that coming from the correction 
to the orbital response. Taking the difference. 



01 



7 



ep 



1/2 



(l + ecosw)2 
X [e(cos u + e) -f 2 sin^ w 
— 2 cos^ v(l + e cost" )] . 

(D17) 

The first term in the square bracket is proportional to the 
Newtonian term and therefore vanishes after averaging. 
The remaining term can be integrated analytically. 



dv' 



siv? v' 



— cos^ v'{l + ecosv') 
(1 -I- ecosw')^ 



sm V cos V 



1 + e cos V 

_(D18) 

which is also zero after integration over one orbit. We 
conclude that the leading order relativistic correction in 
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the perturbation equation predicts the increase in eccen- 
tricity that we observe numerically. This secular change 
comes equally from the first order correction to the ra- 
dial velocity and the first order correction to the orbital 
response to an azimuthal perturbation. The relativistic 
corrections can be thought of as an extra force which 
slows down the effective radial motion and accelerates 
the effective azimuthal motion that enter the drag force. 
The radial drag force is correspondingly reduced, while 
the azimuthal drag force is increased and both drive a 



secular increase in eccentricity. The equality of the two 
parts of the force may reflect some hidden symmetry in 
the equations. The response of the orbit to a pertur- 
bation depends on the velocity at each point along the 
orbit, and we are using that same velocity to prescribe 
the perturbation in this case, which might explain why 
the net contribution from the two terms is equal. How- 
ever, the osculating element equations are not explicit in 
how they depend on the instantaneous velocity, so this is 
only a speculation. 
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